Boundary Value Problem – Axial Temperature Variation of a current-carrying bare wire

(Problem 11.28 from the Numerical Methods for Engineers and Scientists)

This study was developed within the Master’s Degree in Advanced Computing, under the curricular unit Numerical Simulation for Engineering, by Diogo Silva (PG61444) and Tomás Pereira (PG59810).

The objective is to numerically solve the axial temperature distribution along a current-carrying bare wire.

Figure 1 – Axial Temperature Variation of a Current-Carrying Bare Wire

The differential equation that describes the problem is the following:

$$\frac{d^2T}{dx^2} - \frac{4h}{kD}\,(T - T_\infty) - \frac{4 \varepsilon \sigma_{SB}}{kD}\,(T^4 - T_\infty^4) = -\frac{I^2 \rho_e}{k\,(\pi D^2 /4)^2}$$

Boundary Conditions:

• $T(0) = T_\infty$ ($x=0$)
• $\dfrac{dT}{dx}\big|_{x=L/2} = 0$ ($x=\frac{L}{2}$)

Notes:

• We are going to solve the boundary value problem for $0\leq x \leq \frac{L}{2}$.
• The solution is obtained with the help from a solver BVP from the python library scipy.integrate.solve_bvp

1. Import necessary libraries

Firstly we need to import necessary libraries to this notebook.

• numpy is needed for math operations such as numeric arrays, vectorized evaluation of the solution and parameters.
• scipy for the core boundary value problem solver bvp.
• matplotlib for plotting the results.

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

2. Physical Parameters

Below are the physical constants and parameters that define the problem.

These describe the thermal, electrical and geometric characteristics of the current-carrying wire.

# -----------------------------
# Thermal and electrical properties
k = 72.0          # Thermal conductivity [W/(mK)]
h = 2000.0        # Convective heat transfer coefficient [W/(m^2K)]
eps = 0.1         # Surface emissivity
sigma = 5.67e-8   # Stefan–Boltzmann constant [W/(m^2K^4)]
rho_e = 32e-8     # Electrical resistivity [omega m]

# -----------------------------
# Electrical and ambient conditions
I = 2.0           # Electric current [A]
T_inf = 300.0     # Ambient temperature [K]


# -----------------------------
# Geometric parameters
D = 7.62e-5       # Wire diameter [m]
L = 4.0e-3        # Total wire length [m]


A = np.pi * D**2 / 4.0   # [m^2]

3. Auxiliary Constants

To simplify the differential equation, we define some auxiliary constants below. These parameters remain constant trhoughout the domain and allow us to write the differential equation in a more compact form.

C1 = 4*h / (k*D)
C2 = 4*eps*sigma / (k*D)
C3 = I**2 * rho_e / (k*A**2)

print(f"C1 = {C1:.3e}, C2 = {C2:.3e}, C3 = {C3:.3e}")

C1 = 1.458e+06, C2 = 4.134e-06, C3 = 8.548e+08

4. Define the Differential Equation and Boundary Conditions

Second-Order ODE to First-Order System

We started from a second-order differential equation describing the steady-state temperature variation along the wire. $$ \frac{d^2T}{dx^2} - \frac{4h}{kD}(T - T_\infty) - \frac{4\varepsilon\sigma_{SB}}{kD}(T^4 - T_\infty^4) = -\frac{I^2 \rho_e}{k\left(\frac{\pi D^2}{4}\right)^2}. $$

To make it more suitable for numerical solution using solve_bvp, we converted this equation into a system of first-order equations by introducing two new variables:

$$ y_1(x) = T(x), \quad y_2(x) = \frac{dT}{dx}. $$

Substituting these definitions gives:

$$ \begin{cases} y_1' = y_2, \\[4pt] y_2' = C_1 (y_1 - T_\infty) + C_2 (y_1^4 - T_\infty^4) - C_3, \end{cases} $$

where the constants (C_1), (C_2), and (C_3) are defined as above.

Boundary Conditions

Because the temperature distribution is symmetric about the midpoint of the wire, we only need to solve for half the wire length, from (x=0) to (x=L/2).

$$ \begin{cases} y_1(0) = T_\infty = 300\text{ K}, \\[4pt] y_2(L/2) = 0. \end{cases} $$

• The Dirichlet boundary condition fixes the temperature at the wire surface ($x=0$).
• The Neumann boundary condition imposes zero temperature gradient at the midpoint ($x=\frac{L}{2}$), representing symmetry and no heat flux through that plane.

Figure 2 – EDO Resolution

#Function defining the ODE system
# y[0] = T
# y[1] = dT/dx

def fun(x, y):
    T = y[0]
    dTdx = y[1]
    d2Tdx2 = C1*(T - T_inf) + C2*(T**4 - T_inf**4) - C3
    return np.vstack((dTdx, d2Tdx2))

# Boundary conditions
def bc(ya, yb):
    # ya -> values at x=0; yb -> values at x=L/2
    return np.array([
        ya[0] - T_inf,   # T(0) = T_inf
        yb[1]            # dT/dx(L/2) = 0
    ])

5. Initial Guess

The boundary value solver requires an initial estimate for both the temperature $T(x)$ and its derivative $\frac{dT}{dx}$ along the domain.

As a physically reasonable starting point, we assume:

• The wire is initially at the ambient temperature $(T_\infty = 300\,\mathrm{K})$.
• There is no temperature gradient (flat profile).

This provides a simple and stable initial condition for the iterative solver: $$ T(x) = T_\infty, \quad \frac{dT}{dx} = 0. $$

#Initial guess
x = np.linspace(0, L/2, 50)
y_init = np.vstack((
    T_inf * np.ones_like(x),     # palpite inicial T(x)=T_inf
    np.zeros_like(x)             # dT/dx = 0
))

6. BVP Solver

Similar to the MATLAB bvp4c function, we use the solve_bvp function from the scipy.integrate library to solve the boundary value problem numerically.

sol = solve_bvp(fun, bc, x, y_init, tol=1e-6, max_nodes=10000)

if sol.success:
    print("Solution converged.")
else:
    print("Solution did not converge.")

Solution converged.

7. Plot the results

Finally, we plot the temperature distribution along the wire length to visualize the results.

x_plot = np.linspace(0, L/2, 200)
T_plot = sol.sol(x_plot)[0]

# Max temperature at the center (x=L/2)
T_max = T_plot[-1]

# Tentar usar Plotly para interatividade; se não estiver disponível, usar Matplotlib como fallback
try:
    import plotly.graph_objects as go

    hover = 'x: %{x:.2f} mm<br>T: %{y:.2f} K'
    fig = go.Figure()
    fig.add_trace(go.Scatter(
        x=x_plot*1e3,
        y=T_plot,
        mode='lines',
        name='T(x)',
        line=dict(color='firebrick', width=3),
        hovertemplate=hover
    ))

    fig.add_trace(go.Scatter(
        x=[x_plot[-1]*1e3],
        y=[T_max],
        mode='markers+text',
        name='T_max',
        marker=dict(size=8, color='royalblue'),
        text=[f'T_max = {T_max:.2f} K'],
        textposition='top right',
        hovertemplate='x: %{x:.2f} mm<br>T_max: %{y:.2f} K'
    ))

    fig.update_layout(
        title='Distribuição de Temperatura ao longo do fio (0 ≤ x ≤ L/2)',
        xaxis_title='x [mm]',
        yaxis_title='Temperatura T [K]',
        template='plotly_white',
        width=800,
        height=500
    )
    fig.update_xaxes(showgrid=True)
    fig.update_yaxes(showgrid=True)
    fig.show()

except Exception as e:
    # Fallback para Matplotlib
    plt.figure(figsize=(8,5))
    plt.plot(x_plot*1e3, T_plot, 'r', lw=2)
    plt.xlabel('x [mm]')
    plt.ylabel('Temperatura T [K]')
    plt.title('Distribuição de Temperatura ao longo do fio (0 ≤ x ≤ L/2)')
    plt.grid(True)

    # Anotar a temperatura máxima
    plt.scatter([x_plot[-1]*1e3], [T_max], color='blue')
    plt.annotate(f'T_max = {T_max:.2f} K', xy=(x_plot[-1]*1e3, T_max),
                 xytext=(x_plot[-1]*1e3*0.9, T_max + 2),
                 arrowprops=dict(arrowstyle='->', color='black'))
    plt.show()

print(f"T_max no centro (x=L/2): {T_max:.2f} K")

T_max no centro (x=L/2): 781.60 K

8. Conclusion

The steady‑state axial temperature distribution of the current‑carrying bare wire was obtained by solving the boundary value problem on half the length $(0 \leq x \leq \frac{L}{2})$ using symmetry $(T(0)=T_\infty$ and $\frac{dT}{dx}(L/2)=0)$.

The solution rises monotonically from the cooled end to a maximum at the center, where the axial heat flux vanishes. For the given parameters the computed centerline maximum temperature is $T_{max} = 781.60 \, K$.

Reformulating the second‑order ODE as a first‑order system enabled robust convergence with solve_bvp, yielding a smooth temperature field suitable for further parametric or interactive exploration.